        subroutine CC(E_ele,E_nu,Num_ele,NORB,Eorb,ERISPIN,E_cc)
        implicit none
        integer :: kcy
        integer :: ii,ia,ij,ib !spin orb
        integer :: iispa,iaspa,ijspa,ibspa
        integer :: Num_ele,NORB,nspinocc,nspinvir,nspin
        real*8  :: E_ele,E_nu
        real*8  :: Eorb(NORB),Hcore_mo(NORB,NORB),fockspin(2*NORB,2*NORB)
        real*8  :: erispin(2*NORB,2*NORB,2*NORB,2*NORB)
        real*8,allocatable :: tia(:,:),tijab(:,:,:,:),Dia(:,:),Dijab(:,:,:,:)
        real*8,allocatable :: Fvv(:,:),Foo(:,:),Fov(:,:)
        real*8,allocatable :: W_oooo(:,:,:,:),W_vvvv(:,:,:,:),W_ovvo(:,:,:,:)
        real*8,allocatable :: tauijab1(:,:,:,:),tauijab2(:,:,:,:)
        real*8,allocatable :: tianew(:,:),tijabnew(:,:,:,:)
        integer :: i,j,a,b,maxkcy
        real*8  :: Ethreshold,rmsthreshold,conver_Ecc,conver_rmscc
        real*8  :: E_mp2tem,tijabnow,x,E_cc,rms,E_ccold
!!! set diis variable
        integer,Parameter :: ndiis_cc=7
        integer :: inv_DIIS,kq,kp
        real*8,allocatable  :: t1(:,:,:),errmat1(:,:,:)
        real*8,allocatable  :: t2(:,:,:,:,:),errmat2(:,:,:,:,:)
        real*8  :: cof1_cc(ndiis_cc+1),cof2_cc(ndiis_cc+1)
        real*8  :: Bdiis1(ndiis_cc+1,ndiis_cc+1),Bdiis2(ndiis_cc+1,ndiis_cc+1)
        integer ierr,lwork
        real*8,allocatable :: work(:)
        real*8,allocatable :: ipiv(:)
!!!invoke CCSD(T)
        integer :: inv_CCSDT
!!! end
!!! set parameter
        inv_CCSDT=1 !!!=1:calculate CCSD(T) ;/=1:don't invoke calculate
        inv_DIIS=1  !!!=1:invoke DIIS;/=1:don't invoke DIIS
        maxkcy=500
!!!end
        write(21,*)'==========================================='
        write(21,*)'*** Start CCSD(T) procedure ***'
        write(21,'(35X,40A)')'*******************'
        write(21,'(35X,40A)')'* CCSD Iterations *'
        write(21,'(35X,40A)')'*******************'

        if(inv_DIIS==1)then
           write(21,*)'==============='
           write(21,*)'Invoke DIIS_CC'
           write(21,*)'==============='
        endif
        write(21,*)''
!!! Set the maximum number of cycles         
!!!
        nspin=2*NORB
        nspinocc=Num_ele
        nspinvir=nspin-nspinocc
        conver_Ecc=8.0d0
        conver_rmscc=8.0d0
        Ethreshold=10d0**(-conver_Ecc)
        rmsthreshold=10d0**(-conver_rmscc)
        allocate(tia(nspinocc,nspinvir),tijab(nspinocc,nspinocc,nspinvir,nspinvir),&
                tianew(nspinocc,nspinvir),tijabnew(nspinocc,nspinocc,nspinvir,nspinvir),&
                 Dia(nspinocc,nspinvir),Dijab(nspinocc,nspinocc,nspinvir,nspinvir))
        allocate(Fvv(nspinvir,nspinvir),Foo(nspinocc,nspinocc),Fov(nspinocc,nspinvir))
        allocate(W_oooo(nspinocc,nspinocc,nspinocc,nspinocc),&
                 W_vvvv(nspinvir,nspinvir,nspinvir,nspinvir),&
                 W_ovvo(nspinocc,nspinvir,nspinvir,nspinocc))
        allocate(tauijab1(nspinocc,nspinocc,nspinvir,nspinvir),&
                  tauijab2(nspinocc,nspinocc,nspinvir,nspinvir))
        
!!! invoke DIIS
        allocate(t1(nspinocc,nspinvir,ndiis_cc),errmat1(nspinocc,nspinvir,ndiis_cc))
        allocate(t2(nspinocc,nspinocc,nspinvir,nspinvir,ndiis_cc))
        allocate(errmat2(nspinocc,nspinocc,nspinvir,nspinvir,ndiis_cc))
!!!end allocate DIIS
        tia=0.0d0;tijab=0.0d0
        Dia=0.0d0;Dijab=0.0d0        
        !!!! set tia and tijab. Compute MP2.
        E_mp2tem=0d0;Dijab=0.0d0;tijabnow=0.0d0
        do ii=1,nspinocc
            iispa=ceiling(ii/2d0)
            do ia=nspinocc+1,nspin
                iaspa=ceiling(ia/2d0)
                Dia(ii,ia-nspinocc)= Eorb(iispa)-Eorb(iaspa)
                do ij=1,nspinocc
                ijspa=ceiling(ij/2d0)
                    do ib=nspinocc+1,nspin
                    ibspa=ceiling(ib/2d0)
                    Dijab(ii,ij,ia-nspinocc,ib-nspinocc)=Eorb(iispa)+Eorb(ijspa)-Eorb(iaspa)-Eorb(ibspa)
                    tijabnow=ERISPIN(ii,ij,ia,ib)/Dijab(ii,ij,ia-nspinocc,ib-nspinocc)
                    E_mp2tem=E_mp2tem+ERISPIN(ii,ij,ia,ib)*tijabnow
                    tijab(ii,ij,ia-nspinocc,ib-nspinocc)=tijabnow
                    end do
                end do
            end do
        end do
        E_mp2tem=0.25d0*E_mp2tem
        !!!Initial energy E_CC=E_mp2

        E_cc=0.0d0
        do ii=1,nspinocc
           do ia=1,nspinvir
              do ij=1,nspinocc
                 do ib=1,nspinvir
                    !Build the effective singles and doubles operators.
                    x=(tia(ii,ia)*tia(ij,ib)-tia(ii,ib)*tia(ij,ia))
                    tauijab1(ii,ij,ia,ib)=tijab(ii,ij,ia,ib)+0.5d0*x
                    tauijab2(ii,ij,ia,ib)=tauijab1(ii,ij,ia,ib)+0.5d0*x
                    E_cc=E_cc + 0.25d0*tijab(ii,ij,ia,ib)*ERISPIN(ii,ij,ia+nspinocc,ib+nspinocc)
                 enddo
              enddo
           enddo
        enddo
        !write(*,*)E_cc

!!!   for hartree-fock reference       
        fockspin=0.0d0
        do ii=1,nspin
            iispa=ceiling(ii/2d0)
            fockspin(ii,ii)=Eorb(iispa)
        end do

!!!Reference: Stanton et al., J. Chem. Phys. 94, 4334 (1990) (primary)
!!!Reference: Hirata et al., J. Chem. Phys. 120, 2581 (2004)
      write(21,'(2X,A5,3X,6(A10,13X))')'Cycle','E_cc','E_ccold','delt E','rms'
!!!set F_me=F_me +tnf*<mn|ef>  (F_me eqval F_ov)
        do kcy=1,maxkcy
          
           !call F_ov(NORB,nspinocc,nspinvir,tia,fockspin,ERISPIN,Fov)
!!!set F_ae(F_vv)
!F_ae=(1-delta(a,e))*f(a,e)-1/2f(m,e)*tma+tmf*<ma|fe>-1/2&mnaf<mn|ef>
           call F_vv(NORB,nspinocc,nspinvir,tia,fockspin,tauijab1,ERISPIN,Fvv)
!!!set F_mi(F_oo)
           call F_oo(NORB,nspinocc,nspinvir,tia,fockspin,tauijab1,ERISPIN,Foo)
           call F_ov(NORB,nspinocc,nspinvir,tia,fockspin,ERISPIN,Fov)
!!!set W_oooo
           call WW_oooo(NORB,nspinocc,nspinvir,tia,tauijab2,ERISPIN,W_oooo)
!!!set W_vvvv
           call WW_vvvv(NORB,nspinocc,nspinvir,tia,tauijab2,ERISPIN,W_vvvv)
!!!set W_ovvo
           call WW_ovvo(NORB,nspinocc,nspinvir,tia,tijab,ERISPIN,W_ovvo)
!!! update the new tianew
           call update_tia(NORB,nspinocc,nspinvir,tia,tijab,fockspin,Dia,&
                   Fvv,Foo,Fov,ERISPIN,tianew)
!!! update the new tijabnew
           call update_tijab(NORB,nspinocc,nspinvir,tia,tijab,fockspin,Dijab,&
                   Fvv,Foo,Fov,W_oooo,W_vvvv,W_ovvo,tauijab2,ERISPIN,tijabnew)
!!! call diis_cc
        if(inv_DIIS==1)then
                t1(:,:,mod(kcy,ndiis_cc)+1)=tianew
                t2(:,:,:,:,mod(kcy,ndiis_cc)+1)=tijabnew
                errmat1(:,:,mod(kcy,ndiis_cc)+1)=tianew-tia
                errmat2(:,:,:,:,mod(kcy,ndiis_cc)+1)=tijabnew-tijab
         if(kcy>ndiis_cc)then
                Bdiis1=0.0d0;Bdiis2=0.0d0
                do kp=1,ndiis_cc
                   do kq=1,kp
                      Bdiis1(kp,kq)=sum(errmat1(:,:,kp)*errmat1(:,:,kq))
                      Bdiis2(kp,kq)=sum(errmat2(:,:,:,:,kp)*errmat2(:,:,:,:,kq)) 
                   enddo
                enddo
                Bdiis1(ndiis_cc+1,:)=-1.0d0
                Bdiis1(ndiis_cc+1,ndiis_cc+1)=0.0d0
                Bdiis2(ndiis_cc+1,:)=-1.0d0
                Bdiis2(ndiis_cc+1,ndiis_cc+1)=0.0d0
                cof1_cc=0.0d0;cof1_cc(ndiis_cc+1)=-1.0d0
                cof2_cc=0.0d0;cof2_cc(ndiis_cc+1)=-1.0d0   
                !solve New matrix of cof1_cc
                lwork=-1
                allocate(ipiv(size(cof1_cc)))
                do ii = 1,2
                        !lwork=-1
                        allocate(work(abs(lwork)))
                        call dsysv('L', ndiis_cc+1, 1, Bdiis1, ndiis_cc+1,ipiv, cof1_cc,ndiis_cc+1, work, lwork, ierr)
                        lwork = nint(work(1))
                        deallocate(work)
                end do
                if(ierr/=0)stop 'Solve cof1_cc failed!!!'
                deallocate(ipiv)
                !solve New matrix of cof2_cc
                lwork=-1
                allocate(ipiv(size(cof2_cc)))
                do ii = 1,2
                        !lwork=-1
                        allocate(work(abs(lwork)))
                        call dsysv('L', ndiis_cc+1, 1, Bdiis2,ndiis_cc+1,ipiv, cof2_cc,ndiis_cc+1, work, lwork, ierr)
                        lwork = nint(work(1))
                        deallocate(work)
                end do
                if(ierr/=0)stop 'Solve cof2_cc failed!!!'
                deallocate(ipiv)
         
!!!set  ew tia and tijab
                tianew=0.0d0
                tijabnew=0.0d0
                do i=1,ndiis_cc
                        tianew=tianew+cof1_cc(i)*t1(:,:,i)
                        tijabnew=tijabnew+cof2_cc(i)*t2(:,:,:,:,i)
                enddo 
         endif
        endif  !!! inv_DIIS
!!! end DIIS 


!!! update Energy
           !E_ccold=E_cc
           E_cc=0.0d0
          
           do i=1,nspinocc
              do a=1,nspinvir
                 E_cc=E_cc+fockspin(i,a+nspinocc)*tianew(i,a)
                 do j=1,nspinocc
                    do b=1,nspinvir
                       E_cc=E_cc+0.25d0*tijabnew(i,j,a,b)*ERISPIN(i,j,a+nspinocc,b+nspinocc)
                       E_cc=E_cc+0.50d0*tianew(i,a)*tianew(j,b)*ERISPIN(i,j,a+nspinocc,b+nspinocc)
                    enddo
                 enddo
              enddo
           enddo         
!!! calculate Rms
           rms=sum(sqrt((tia-tianew)**2))
           write(21,'(I4,5(F23.16))')kcy,E_cc,E_ccold,E_cc-E_ccold,rms
           if(abs(E_cc-E_ccold)<=Ethreshold .and. abs(rms)<=rmsthreshold)then 
               write(21,*)'==========================================='
               write(21,*)'Congratulations, CCSD has converged!'
               write(21,*)'Number of iterations:',kcy
               write(21,*)'The CC convergence threshold of E:',Ethreshold,'A.U.'
               write(21,*)'E_CCSD:            ',E_cc
               exit
             else if(kcy==maxkcy)then
                write(21,*) 'Waring: Bad convergence of CC!'
                write(21,*) 'Waring: Bad convergence of CC!'
                write(21,*)'==========================================='    
           endif
           E_ccold=E_cc
           tia=0.0d0;tijab=0.0d0
           tia=tianew
           tijab=tijabnew
           do ii=1,nspinocc
               do ia=1,nspinvir
                       !write(25,*)ii,ia,tia(ii,ia)
                   do ij=1,nspinocc
                       do ib=1,nspinvir
                        x=(tia(ii,ia)*tia(ij,ib)-tia(ii,ib)*tia(ij,ia))
                        tauijab1(ii,ij,ia,ib)=tijab(ii,ij,ia,ib)+0.5d0*x
                        tauijab2(ii,ij,ia,ib)=tauijab1(ii,ij,ia,ib)+0.5d0*x
                        !write(25,*)ii,ij,ia,ib,tauijab1(ii,ij,ia,ib)
                       end do
                   end do
               end do
           end do
        end do  !!!kcy
        write(21,*)'Energr total(CCSD):',(E_nu+E_ele)+E_cc
        write(21,*)'*** The CCSD procedure Done **'
        write(21,*)'==========================================='
!!!CCSD(T) 
        if(inv_CCSDT==1)then
        call CCSD_T(E_cc,E_ele,E_nu,NORB,nspinocc,nspinvir,tia,tijab,Eorb,ERISPIN)
        endif
!!!  deallocate array
        deallocate(tia,tijab,tianew,tijabnew,Dia,Dijab)
        deallocate(Fvv,Foo,Fov)
        deallocate(W_oooo,W_vvvv,W_ovvo)
        deallocate(tauijab1,tauijab2)

!!! DIIS array
        deallocate(t2,errmat2)

        endsubroutine CC
        
        subroutine F_ov(NORB,nspinocc,nspinvir,tia,fockspin,ERISPIN,F_me)
        implicit none
        integer :: NORB
        integer :: m,n,e,f  !occ:m,n; vir:e,f
        integer :: nspinocc,nspinvir
        real*8  :: fockspin(2*NORB,2*NORB)
        real*8  :: F_me(nspinocc,nspinvir),tia(nspinocc,nspinvir)
        real*8  :: ERISPIN(2*NORB,2*NORB,2*NORB,2*NORB)
        real*8  :: tem 
        F_me=0.0d0
        tem =0.0d0
        do m=1,nspinocc
           do e=1,nspinvir
              F_me(m,e)=fockspin(m,e+nspinocc)
              tem=0.0d0
              do n=1,nspinocc
                 do f=1,nspinvir
                    tem=tem+tia(n,f)*ERISPIN(m,n,e+nspinocc,f+nspinocc)
                 enddo
              enddo
              F_me(m,e)=F_me(m,e)+tem
           enddo
        enddo
        endsubroutine F_ov

        subroutine F_vv(NORB,nspinocc,nspinvir,tia,fockspin,tauijab1,ERISPIN,F_ae)
        implicit none
        integer :: NORB
        integer :: m,n,e,f,a
        integer :: nspinocc,nspinvir
        real*8  :: fockspin(2*NORB,2*NORB)
        real*8  :: F_ae(nspinvir,nspinvir),tia(nspinocc,nspinvir)
        real*8  :: ERISPIN(2*NORB,2*NORB,2*NORB,2*NORB)
        real*8  :: tauijab1(nspinocc,nspinocc,nspinvir,nspinvir)
        real*8  :: krodelta
        real*8  :: tem
        F_ae=0.0d0
        tem =0.0d0
        do a=1,nspinvir
           do e=1,nspinvir
              !!!first
              F_ae(a,e)=(1.0d0-krodelta(a,e))*fockspin(a+nspinocc,e+nspinocc)
              !!!second 
              tem=0.0d0
              do m=1,nspinocc
                tem=tem + fockspin(m,e+nspinocc)*tia(m,a)
              enddo
              F_ae(a,e)=F_ae(a,e)-0.50d0*tem 
              !!!third 
              tem=0.0d0
              do m=1,nspinocc
                 do f=1,nspinvir
                    tem=tem+tia(m,f)*ERISPIN(m,a+nspinocc,f+nspinocc,e+nspinocc)
                 enddo
              enddo
              F_ae(a,e)=F_ae(a,e)+tem
              !!!fourth
              tem=0.0d0
              do m=1,nspinocc
                 do n=1,nspinocc
                    do f=1,nspinvir
                    tem=tem+tauijab1(m,n,a,f)*ERISPIN(m,n,e+nspinocc,f+nspinocc)    
                    enddo
                 enddo
              enddo
              F_ae(a,e)=F_ae(a,e)-0.50d0*tem
           enddo 
        enddo
        endsubroutine F_vv       
                 

        subroutine F_oo(NORB,nspinocc,nspinvir,tia,fockspin,tauijab1,ERISPIN,F_mi)
        implicit none
        integer :: NORB
        integer :: m,n,e,f,i
        integer :: nspinocc,nspinvir
        real*8  :: fockspin(2*NORB,2*NORB)
        real*8  :: F_mi(nspinocc,nspinocc),tia(nspinocc,nspinvir)
        real*8  :: ERISPIN(2*NORB,2*NORB,2*NORB,2*NORB)
        real*8  :: tauijab1(nspinocc,nspinocc,nspinvir,nspinvir)
        real*8  :: krodelta
        real*8  :: tem
        F_mi=0.0d0
        tem=0.0d0
        do m=1,nspinocc
           do i=1,nspinocc
              !!!!first
              F_mi(m,i)=(1.0d0-krodelta(m,i))*fockspin(m,i)
              !!!second
              tem=0.0d0
              do e=1,nspinvir
                 tem=tem+tia(i,e)*fockspin(m,e+nspinocc)
              enddo
              F_mi(m,i)=F_mi(m,i)+0.50d0*tem
              !!!third
              tem=0.0d0
              do e=1,nspinvir
                 do n=1,nspinocc
                  tem=tem+tia(n,e)*ERISPIN(m,n,i,e+nspinocc)
                 enddo
              enddo
              F_mi(m,i)=F_mi(m,i)+ tem
              !!!fourth
              tem=0.0d0
              do n=1,nspinocc
                 do e=1,nspinvir
                    do f=1,nspinvir
                     tem=tem+tauijab1(i,n,e,f)*ERISPIN(m,n,e+nspinocc,f+nspinocc)
                    enddo       
                 enddo
              enddo
              F_mi(m,i)=F_mi(m,i)+0.50d0*tem
           enddo
        enddo 

        endsubroutine F_oo

        subroutine WW_oooo(NORB,nspinocc,nspinvir,tia,tauijab2,ERISPIN,Woooo)
        implicit none
        integer :: m,n,i,j,e,f  !!occ :m,n,i,j ; vir :i,j
        integer :: NORB
        integer :: nspinocc,nspinvir
        real*8  :: tia(nspinocc,nspinvir)
        real*8  :: ERISPIN(2*NORB,2*NORB,2*NORB,2*NORB)
        real*8  :: tauijab2(nspinocc,nspinocc,nspinvir,nspinvir)
        real*8  :: Woooo(nspinocc,nspinocc,nspinocc,nspinocc)
        real*8  :: krodelta
        real*8  :: tem
        Woooo=0.0d0
        tem=0.0d0
        do m=1,nspinocc
           do n=1,nspinocc
              do i=1,nspinocc
                 do j=1,nspinocc
                 !!!first
                 Woooo(m,n,i,j)=ERISPIN(m,n,i,j)
                 tem=0.0d0
                 !!!second
                 do e=1,nspinvir
                    tem=tem+tia(j,e)*ERISPIN(m,n,i,e+nspinocc)-tia(i,e)*ERISPIN(m,n,j,e+nspinocc)
                 enddo
                 Woooo(m,n,i,j)=Woooo(m,n,i,j)+tem
                 !!!third
                 tem=0.0d0
                 do e=1,nspinvir
                    do f=1,nspinvir
                      tem=tem+tauijab2(i,j,e,f)*ERISPIN(m,n,e+nspinocc,f+nspinocc)
                    enddo
                 enddo
                 Woooo(m,n,i,j)=Woooo(m,n,i,j)+0.25d0*tem
                 enddo
              enddo
           enddo
        enddo
        endsubroutine WW_oooo

        subroutine WW_vvvv(NORB,nspinocc,nspinvir,tia,tauijab2,ERISPIN,Wvvvv)
        implicit none
        integer :: a,b,e,f,m,n  !!occ :m,n; vir :a,b,e,f
        integer :: NORB
        integer :: nspinocc,nspinvir
        real*8  :: tia(nspinocc,nspinvir)
        real*8  :: ERISPIN(2*NORB,2*NORB,2*NORB,2*NORB)
        real*8  :: tauijab2(nspinocc,nspinocc,nspinvir,nspinvir)
        real*8  :: Wvvvv(nspinvir,nspinvir,nspinvir,nspinvir)
        real*8  :: krodelta
        real*8  :: tem
        Wvvvv=0.0d0
        do a=1,nspinvir
           do b=1,nspinvir
              do e=1,nspinvir
                 do f=1,nspinvir
                 !!!first
                 Wvvvv(a,b,e,f)=ERISPIN(a+nspinocc,b+nspinocc,e+nspinocc,f+nspinocc)
                 !!!second
                 tem=0.0d0
                 do m=1,nspinocc
                    tem=tem+tia(m,b)*ERISPIN(a+nspinocc,m,e+nspinocc,f+nspinocc)&
                        -tia(m,a)*ERISPIN(b+nspinocc,m,e+nspinocc,f+nspinocc)
                 enddo
                Wvvvv(a,b,e,f)=Wvvvv(a,b,e,f)-tem
                !!!third
                tem=0.0d0
                do m=1,nspinocc
                   do n=1,nspinocc
                      tem=tem+tauijab2(m,n,a,b)*ERISPIN(m,n,e+nspinocc,f+nspinocc)
                   enddo
                enddo
                Wvvvv(a,b,e,f)=Wvvvv(a,b,e,f)+0.25d0*tem
                enddo
             enddo
           enddo
        enddo
        endsubroutine WW_vvvv   

        subroutine WW_ovvo(NORB,nspinocc,nspinvir,tia,tijab,ERISPIN,Wovvo)
        implicit none
        integer :: m,b,e,j,n,f  !!occ :m,j,n; vir :b,e,f
        integer :: NORB
        integer :: nspinocc,nspinvir
        real*8  :: tia(nspinocc,nspinvir)
        real*8  :: ERISPIN(2*NORB,2*NORB,2*NORB,2*NORB)
        real*8  :: tijab(nspinocc,nspinocc,nspinvir,nspinvir)
        real*8  :: Wovvo(nspinocc,nspinvir,nspinvir,nspinocc)
        real*8  :: krodelta
        real*8  :: tem
        Wovvo=0.0d0
        do m=1,nspinocc
           do b=1,nspinvir  
              do e=1,nspinvir
                 do j=1,nspinocc  
                  !!!first
                  Wovvo(m,b,e,j)=ERISPIN(m,b+nspinocc,e+nspinocc,j)
                  !!!second
                  tem=0.0d0
                  do f=1,nspinvir
                     tem=tem+tia(j,f)*ERISPIN(m,b+nspinocc,e+nspinocc,f+nspinocc)
                  enddo
                  Wovvo(m,b,e,j)=Wovvo(m,b,e,j)+tem
                  !!!third
                  tem=0.0d0
                  do n=1,nspinocc
                     tem=tem+tia(n,b)*ERISPIN(m,n,e+nspinocc,j)
                  enddo
                  Wovvo(m,b,e,j)=Wovvo(m,b,e,j)-tem
                  !!!fourth
                  tem=0.0d0
                  do n=1,nspinocc
                     do f=1,nspinvir
                        tem=tem+(0.5d0*tijab(j,n,f,b)+tia(j,f)*tia(n,b))*ERISPIN(m,n,e+nspinocc,f+nspinocc)
                     enddo
                  enddo
                  Wovvo(m,b,e,j)=Wovvo(m,b,e,j)-tem
                  enddo
               enddo
            enddo
        enddo
        endsubroutine WW_ovvo

        subroutine update_tia(NORB,nspinocc,nspinvir,tia,tijab,fockspin,Dia,Fae,Fmi,Fme,ERISPIN,tianew)
        implicit none
        integer :: NORB
        integer :: i,a,e,m,n,f
        integer :: nspinocc,nspinvir
        real*8  :: Dia(nspinocc,nspinvir)
        real*8  :: fockspin(2*NORB,2*NORB)
        real*8  :: tia(nspinocc,nspinvir),tijab(nspinocc,nspinocc,nspinvir,nspinvir)
        real*8  :: Fme(nspinocc,nspinvir),Fmi(nspinocc,nspinocc)
        real*8  :: Fae(nspinvir,nspinvir)
        real*8  :: ERISPIN(2*NORB,2*NORB,2*NORB,2*NORB)
        real*8  :: tianew(nspinocc,nspinvir)
        real*8  :: tem,tianow
        tianow=0.0d0
        tianew=0.0d0
        do i=1,nspinocc
           do a=1,nspinvir
              !!!first
              tianow=0.0d0
              tianow=fockspin(i,a+nspinocc)
              !!!second
              tem=0.0d0
              do e=1,nspinvir
                 tem=tem+tia(i,e)*Fae(a,e)
              enddo
              tianow=tianow+tem
              !!!third
              tem=0.0d0
              do m=1,nspinocc
                 tem=tem+tia(m,a)*Fmi(m,i)
              enddo
              tianow=tianow-tem
              !!!fourth
              tem=0.0d0
              do m=1,nspinocc
                 do e=1,nspinvir
                    tem=tem+tijab(i,m,a,e)*Fme(m,e)
                 enddo
              enddo
              tianow=tianow+tem
              !!!fifth
              tem=0.0d0
              do n=1,nspinocc
                 do f=1,nspinvir
                    tem=tem+tia(n,f)*ERISPIN(n,a+nspinocc,i,f+nspinocc)
                 enddo
              enddo
              tianow=tianow-tem
               !write(25,*)i,a,tianow
              !!! sixth
              tem=0.0d0
              do m=1,nspinocc
                 do e=1,nspinvir
                    do f=1,nspinvir
                       tem=tem+tijab(i,m,e,f)*ERISPIN(m,a+nspinocc,e+nspinocc,f+nspinocc)
                    enddo
                 enddo
              enddo
              tianow=tianow-0.5d0*tem
              !write(25,*)i,a,tianow
              !!!seventh
              tem=0.0d0
              do m=1,nspinocc
                 do e=1,nspinvir
                    do n=1,nspinocc
                       tem=tem+tijab(m,n,a,e)*ERISPIN(n,m,e+nspinocc,i)
                       !write(25,*)m,n,a,e,tijab(m,n,a,e)
                       !write(25,*)n,e+nspinocc,m,i,ERISPIN(n,m,e+nspinocc,i)
                    enddo
                 enddo
              enddo
              tianow=tianow-0.5d0*tem
              !write(25,*)i,a,tianow
              tianew(i,a)=tianow/Dia(i,a)
           enddo
        enddo
        endsubroutine update_tia
              
        subroutine update_tijab(NORB,nspinocc,nspinvir,tia,tijab,fockspin,Dijab,&
                   Fae,Fmi,Fme,W_oooo,W_vvvv,W_ovvo,tauijab2,ERISPIN,tijabnew) 
        implicit none
        integer :: NORB
        integer :: i,j,a,b,e,m,n,f
        integer :: nspinocc,nspinvir
        real*8  :: Dijab(nspinocc,nspinocc,nspinvir,nspinvir)
        real*8  :: fockspin(2*NORB,2*NORB)
        real*8  :: tia(nspinocc,nspinvir),tijab(nspinocc,nspinocc,nspinvir,nspinvir)
        real*8  :: Fme(nspinocc,nspinvir),Fmi(nspinocc,nspinocc)
        real*8  :: Fae(nspinvir,nspinvir)
        real*8  :: ERISPIN(2*NORB,2*NORB,2*NORB,2*NORB)
        real*8  :: tijabnew(nspinocc,nspinocc,nspinvir,nspinvir)
        real*8  :: tauijab2(nspinocc,nspinocc,nspinvir,nspinvir)
        real*8  :: W_oooo(nspinocc,nspinocc,nspinocc,nspinocc)!Wmnij
        real*8  :: W_vvvv(nspinvir,nspinvir,nspinvir,nspinvir)!Wabef
        real*8  :: W_ovvo(nspinocc,nspinvir,nspinvir,nspinocc)!Wmbej
        real*8  :: tem,tem1,tem2,tem3,tem4,tijabnow
        tijabnew=0.0d0
        tijabnow=0.0d0
        do i=1,nspinocc
           do a=1,nspinvir
              do j=1,nspinocc
                 do b=1,nspinvir
                    !!!first
                    tijabnow=ERISPIN(i,j,a+nspinocc,b+nspinocc)
                    !!!second
                    tem=0.0d0
                    do e=1,nspinvir
                       tem1=0.0d0
                       tem2=0.0d0
                       do m=1,nspinocc
                       tem1=tem1+tia(m,b)*Fme(m,e)
                       tem2=tem2+tia(m,a)*Fme(m,e)
                       enddo
                       tem1=tem1*0.5d0
                       tem2=tem2*0.5d0 
                       tem=tem+tijab(i,j,a,e)*(Fae(b,e)-tem1)-tijab(i,j,b,e)*(Fae(a,e)-tem2)
                    enddo
                    tijabnow=tijabnow+tem
                    !!!third
                    tem=0.0d0
                    do m=1,nspinocc
                       tem1=0.0d0
                       tem2=0.0d0
                       do e=1,nspinvir
                          tem1=tem1+tia(j,e)*Fme(m,e)
                          tem2=tem2+tia(i,e)*Fme(m,e)
                       enddo    
                       tem1=tem1*0.5d0
                       tem2=tem2*0.5d0
                       tem=tem+tijab(i,m,a,b)*(Fmi(m,j)+tem1)-tijab(j,m,a,b)*(Fmi(m,i)+tem2)   
                    enddo
                    tijabnow=tijabnow-tem
                    !!!fourth
                    tem=0.0d0
                    do m=1,nspinocc
                       do n=1,nspinocc
                          tem=tem+tauijab2(m,n,a,b)*W_oooo(m,n,i,j)
                        enddo
                    enddo
                    tijabnow=tijabnow+0.5d0*tem
                    !!!fifth
                    tem=0.0d0
                    do e=1,nspinvir
                       do f=1,nspinvir
                          tem=tem+tauijab2(i,j,e,f)*W_vvvv(a,b,e,f)
                       enddo
                    enddo
                    tijabnow=tijabnow+0.5d0*tem
                    !!!sixth
                    tem=0.0d0
                    do m=1,nspinocc
                      !tem1=0.0d0
                      !tem2=0.0d0
                      !tem3=0.0d0
                      !tem4=0.0d0
                       do e=1,nspinvir
                          tem1=tijab(i,m,a,e)*W_ovvo(m,b,e,j)-tia(i,e)*tia(m,a)*&
                               ERISPIN(m,b+nspinocc,e+nspinocc,j)
                          tem2=tijab(i,m,b,e)*W_ovvo(m,a,e,j)-tia(i,e)*tia(m,b)*&
                               ERISPIN(m,a+nspinocc,e+nspinocc,j)
                          tem3=tijab(j,m,a,e)*W_ovvo(m,b,e,i)-tia(j,e)*tia(m,a)*&
                               ERISPIN(m,b+nspinocc,e+nspinocc,i)
                          tem4=tijab(j,m,b,e)*W_ovvo(m,a,e,i)-tia(j,e)*tia(m,b)*&
                               ERISPIN(m,a+nspinocc,e+nspinocc,i)
                          tem=tem+tem1-tem2-tem3+tem4     
                       enddo
                    enddo
                    tijabnow=tijabnow+tem
                    !!!seventh
                    tem=0.0d0
                    do e=1,nspinvir
                      !tem1=0.0d0
                      !tem2=0.0d0
                       tem1=tia(i,e)*ERISPIN(a+nspinocc,b+nspinocc,e+nspinocc,j)
                       tem2=tia(j,e)*ERISPIN(a+nspinocc,b+nspinocc,e+nspinocc,i)
                       tem=tem+tem1-tem2
                    enddo    
                    tijabnow=tijabnow+tem
                    !!!!eighth
                    tem=0.0d0
                    do m=1,nspinocc
                       tem1=0.0d0
                       tem2=0.0d0
                       tem1=tia(m,a)*ERISPIN(m,b+nspinocc,i,j)
                       tem2=tia(m,b)*ERISPIN(m,a+nspinocc,i,j)
                       tem=tem+tem1-tem2
                    enddo
                    tijabnow=tijabnow-tem
                    tijabnew(i,j,a,b)=tijabnow/Dijab(i,j,a,b)
                    !write(25,*)i,j,a,b,tijabnow
                 enddo
              enddo
           enddo
        enddo
        endsubroutine update_tijab  
        subroutine CCSD_T(E_cc,E_ele,E_nu,NORB,nspinocc,nspinvir,tia,tijab,Eorb,ERISPIN)
        implicit none
        integer :: NORB,nspinocc,nspinvir
        integer :: i,j,k,m,a,b,c,e  !!occ:i,j,k,m; vir: a,b,c e
        integer :: ispa,jspa,kspa,mspa,aspa,bspa,cspa,espa
        real*8  :: E_ele,E_nu,E_cc
        real*8  :: Eorb(NORB)
        real*8  :: ERISPIN(2*NORB,2*NORB,2*NORB,2*NORB)
        real*8  :: Eccsdt,tem,tem1,tem2,D3o3v,E_tem,temp
        real*8  :: tia(nspinocc,nspinvir),tijab(nspinocc,nspinocc,nspinvir,nspinvir)
        write(21,*)''
        write(21,*)'*** Start CCSD(T) procedure ***'
        Eccsdt=0.0d0
        !write(*,*)Eorb
        E_tem=0.0d0
        tem2=0.0d0
        do i=1,nspinocc
           ispa=ceiling(i/2.0d0)
           do j=1,nspinocc
              jspa=ceiling(j/2.0d0)
              do k=1,nspinocc
                 kspa=ceiling(k/2.0d0)
                 do a=1,nspinvir
                    aspa=ceiling((a+nspinocc)/2.0d0)
                    do b=1,nspinvir
                       bspa=ceiling((b+nspinocc)/2.0d0)
                       do c=1,nspinvir
                          cspa=ceiling((c+nspinocc)/2.0d0)
                          tem=0.0d0
                          tem1=0.0d0
                          D3o3v=Eorb(ispa)+Eorb(jspa)+Eorb(kspa)-(Eorb(aspa)+Eorb(bspa)+Eorb(cspa))
                          do e=1,nspinvir
                             tem=tem+(tijab(j,k,a,e)*ERISPIN(e+nspinocc,i,b+nspinocc,c+nspinocc)-&
                                tijab(j,k,b,e)*ERISPIN(e+nspinocc,i,a+nspinocc,c+nspinocc)-&
                                tijab(j,k,c,e)*ERISPIN(e+nspinocc,i,b+nspinocc,a+nspinocc))-&
                               (tijab(i,k,a,e)*ERISPIN(e+nspinocc,j,b+nspinocc,c+nspinocc)-&
                                tijab(i,k,b,e)*ERISPIN(e+nspinocc,j,a+nspinocc,c+nspinocc)-&
                                tijab(i,k,c,e)*ERISPIN(e+nspinocc,j,b+nspinocc,a+nspinocc))-&
                               (tijab(j,i,a,e)*ERISPIN(e+nspinocc,k,b+nspinocc,c+nspinocc)-&
                                tijab(j,i,b,e)*ERISPIN(e+nspinocc,k,a+nspinocc,c+nspinocc)-&
                                tijab(j,i,c,e)*ERISPIN(e+nspinocc,k,b+nspinocc,a+nspinocc))
                          enddo!e
                          do m=1,nspinocc
                             tem1=tem1+(tijab(i,m,b,c)*ERISPIN(m,a+nspinocc,j,k)-&
                                        tijab(i,m,a,c)*ERISPIN(m,b+nspinocc,j,k)-&
                                        tijab(i,m,b,a)*ERISPIN(m,c+nspinocc,j,k))-&
                                       (tijab(j,m,b,c)*ERISPIN(m,a+nspinocc,i,k)-&
                                        tijab(j,m,a,c)*ERISPIN(m,b+nspinocc,i,k)-&
                                        tijab(j,m,b,a)*ERISPIN(m,c+nspinocc,i,k))-&
                                       (tijab(k,m,b,c)*ERISPIN(m,a+nspinocc,j,i)-&
                                        tijab(k,m,a,c)*ERISPIN(m,b+nspinocc,j,i)-&
                                        tijab(k,m,b,a)*ERISPIN(m,c+nspinocc,j,i)) 
                          enddo!m
!!!D_ijk^abc*t_ijk^abc(d)
                          tem2=(tia(i,a)*ERISPIN(j,k,b+nspinocc,c+nspinocc)-&
                                tia(i,b)*ERISPIN(j,k,a+nspinocc,c+nspinocc)-&
                                tia(i,c)*ERISPIN(j,k,b+nspinocc,a+nspinocc))-&
                               (tia(j,a)*ERISPIN(i,k,b+nspinocc,c+nspinocc)-&
                                tia(j,b)*ERISPIN(i,k,a+nspinocc,c+nspinocc)-&
                                tia(j,c)*ERISPIN(i,k,b+nspinocc,a+nspinocc))-&
                               (tia(k,a)*ERISPIN(j,i,b+nspinocc,c+nspinocc)-&
                                tia(k,b)*ERISPIN(j,i,a+nspinocc,c+nspinocc)-&
                                tia(k,c)*ERISPIN(j,i,b+nspinocc,a+nspinocc))
                       temp=tem-tem1
                       E_tem=E_tem+temp/D3o3v*(temp+tem2)
                       enddo!c
                    enddo!b
                 enddo!a
              enddo !k
           enddo!j
        enddo!i
        Eccsdt=E_tem/36.0d0
        write(21,'(25X,A30)')'*****  CCSD(T) Summary *****'
        write(21,*) 'CCSD Final results (a.u.)'
        write(21,*)'E(HF_tol)  :',E_ele+E_nu
        write(21,*)'E(CCSD)    :',E_cc
        write(21,*)'E(HF+ECCSD):',E_ele+E_nu+E_cc
        write(21,*)'E(CCSD(T)) :',Eccsdt
        write(21,*)'E(HF+ECCSD+CCSD(T))  :',E_ele+E_nu+E_cc+Eccsdt
        write(21,*)'*** The CCSD(T) procedure Done **'
        endsubroutine CCSD_T
        
